Analyzing the linear frequency modulation using the Wigner-Ville distribution

In [1]:

    

import DSP
using PyPlot


    

In [2]:

    

# include all modules in juwvid
include("../juwvid.jl")


    

    
Out[2]:





juwvid

In [3]:

    

nsample=1024
x,y,iw,ynorm=sampledata.genlinfm(nsample,1.0,0.01);
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(x,y)
PyPlot.savefig("linearFM.png")


    

In [4]:

    

### generating the analytic signal of y
zs=DSP.Util.hilbert(y);


    

In [5]:

    

tfrs=cohenclass.tfrwv(zs);


    

    




Single Wigner Ville
Use fft.

In [6]:

    

a=juwplot.tfrshow(tfrs,x[2]-x[1],x[1],x[end])
#colorbar(a) # uncomment if you show a colorbar
PyPlot.xlabel("time")
PyPlot.ylabel("Frequency")


    

    













    
Out[6]:





PyObject Text(24,0.5,'Frequency')

In [7]:

    

### extracting the instantaneous frequency (IF) from the TFR
indfs=extif.maxif(abs.(tfrs));
dx=x[2]-x[1]
if_ffts=juwutils.index_to_frequency(indfs,NaN,dx,nsample);


    

    




Assuming nft = nsample.

In [8]:

    

### displaying the extracted IF
PyPlot.plot(x,if_ffts,color="C0",lw=3)
PyPlot.plot(x,iw/(2*pi),ls="dashed",color="C3",lw=3)
PyPlot.ylabel("Instantaneous Frequency")
PyPlot.xlabel("Time")


    

    













    
Out[8]:





PyObject Text(0.5,24,'Time')

In [9]:

    

nsample=1024
x,y,iw,ynorm=sampledata.genlinfm(nsample,1.0,0.0001);
PyPlot.plot(x,y)


    

    













    
Out[9]:





1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7fcab5637d30>

In [10]:

    

### z is the analytic signal of y
z=DSP.Util.hilbert(y);


    

In [11]:

    

tfr=cohenclass.tfrwv(z);


    

    




Single Wigner Ville
Use fft.

In [12]:

    

a=juwplot.tfrshow(tfr,x[2]-x[1],x[1],x[end])
PyPlot.xlabel("time")
PyPlot.ylabel("Frequency")
# Because the amplitude is small, you hardly see the modulation by eye.


    

    













    
Out[12]:





PyObject Text(24,0.5,'Frequency')

In [13]:

    

# If you zoom up the TFR, you find that the resolution is poor because of the sampling rate.
a=juwplot.tfrshow(tfr,x[2]-x[1],x[1],x[end],NaN,NaN,15.0)
PyPlot.ylim(0.14,0.18)
PyPlot.xlabel("time")
PyPlot.ylabel("Frequency")


    

    













    
Out[13]:





PyObject Text(24,0.5,'Frequency')

In [14]:

    

### As shown in this figure, the FFT(DFT) cannot provide the sufficient resolution of the frequency.
indf=extif.maxif(abs.(tfr));
dx=x[2]-x[1]
if_fft=juwutils.index_to_frequency(indf,NaN,dx,nsample);
PyPlot.plot(x,if_fft,color="blue",lw=3)
PyPlot.plot(x,iw/(2*pi),ls="dashed",color="red",lw=3)
PyPlot.ylabel("Instantaneous Frequency")
PyPlot.xlabel("Time")
PyPlot.ylim(0.158,0.165)


    

    













    




Assuming nft = nsample.






    
Out[14]:





(0.158, 0.165)

In [15]:

    

### Setting new frequency grid by dividing the frequency indices of f = 90 - 110 by 1024 grids
fin=collect(linspace(90,110,1024));
### Then, we compute the Wigner Ville distribution with the non-uniform FFT.
tfrnufft=cohenclass.tfrwv(z,NaN,NaN,fin,NaN,0);


    

    




Single Wigner Ville
Use nufft.

In [16]:

    

a=juwplot.tfrshow(real(tfrnufft),x[2]-x[1],x[1],x[end],fin[1],fin[end],0.5)
PyPlot.xlabel("time")
PyPlot.ylabel("Frequency")


    

    













    
Out[16]:





PyObject Text(24,0.5,'Frequency')

In [17]:

    

indnf=extif.maxif(abs.(tfrnufft));
dx=x[2]-x[1]
if_nufft=juwutils.index_to_frequency(indnf,fin,dx,nsample);


    

    




Assuming nft = nsample.

In [18]:

    

PyPlot.plot(x,if_fft,color="blue",lw=3)
PyPlot.plot(x,if_nufft,color="green",lw=3)
PyPlot.plot(x,iw/(2*pi),ls="dashed",color="red",lw=3)
PyPlot.ylabel("Instantaneous Frequency")
PyPlot.xlabel("Time")
PyPlot.ylim(0.158,0.165)


    

    













    
Out[18]:





(0.158, 0.165)

In [ ]: